Wmtil 

kSm 


Wmi 


m 

KH&BRrSBESwi 


Hbb 


LIBRARY  OF  THE 

UNIVERSITY  OF  ILLINOIS 

AT  URBANA-CHAMPAIGN 

5/0.84 

no.  S9S-&&0 
cof>2. 


30 


5/c  *r 

r./     'UIUCDCS-R-T3-598 


ytyu^l 


COO-2383-0001 


1 


ASYMPTOTIC  ESTIMATION   OF  ERRORS  AND  DERIVATIVES   FOR   THE  NUMERICAL 
SOLUTION   OF  ORDINARY   DIFFERENTIAL  EQUATIONS 

by 

C.    W.    Gear 
October  1973 


DEPARTMENT  OF  COMPUTER  SCIENCE 
UNIVERSITY  OF  ILLINOIS  AT  URBANA-CHAMPAIGN 


URBANA,  ILLINOIS 


JHE  LIBRARY. OF  THE 
Mftti 


Digitized  by  the  Internet  Archive 
in  2013 


http://archive.org/details/asymptoticestima598gear 


UIUCDCS-R-73-598 


ASYMPTOTIC  ESTIMATION   OF  ERRORS   AND  DERIVATIVES   FOR   THE  NUMERICAL 
SOLUTION   OF   ORDINARY   DIFFERENTIAL  EQUATIONS* 

by 
C.   W.    Gear 


October  1973 


DEPARTMENT  OF   COMPUTER  SCIENCE 

UNIVERSITY   OF  ILLINOIS  AT   URBANA-CHAMPAIGN 

URBANA,    ILLINOIS    6l801 


*     Supported  in  part  by   contract   US  AEC  AT( 11-1)2  383 
PREPRINT   FOR  LIMITED  DISTRIBUTION.      SUBMITTED  TO   IFIPs    "jk. 


ABSTRACT 

Shampine    [h]   points   out   that  the  Milne   device   for  estimating 
errors   may   give   asymptotically  incorrect   results.      This  paper  examines 
general  classes   of  error  estimates    and  shows  how  they   can  he  used  to 
form  asymptotically   correct  estimates   of  the   local  error  in  the 
numerical  solution   of  ordinary  differential  equations,   or  to  estimate 
the  derivatives   of  the  solution  of  the   differential  equation. 


1.      INTRODUCTION 

Shampine    [h]   gives  the   following  example   to   show  that  Milne's 
device   for  estimating  the  error  is  not   asymptotically   correct  in 
some   cases. 
Predictor:      Pn+2   =  yn   +  2h   f(yQ+1)  (D 

Corrector:      Y^2  =  yn+1  +  |(f(7n+1)    +  f(yn+2>]  (2) 

(Asymptotically  it    does   not   matter  whether  we  solve  the  corrector  in  eq. 
(2)    exactly  or  iterate   it    one   or  more  times.)      If  eqs.    (l)    and   (2)    are 
applied  to  y'    =  y,  y(0)   =   1  using  the  exact   initial  conditions    for 
a  fixed  step  size   of  y     =  0,  y     =  e    ,  we    find  that 

p2  -y2  =  -3jh3  +  o(hU)  (3) 

Pm  -  ym  =   -  \  h3  +  0(hU)  m  >_  3  (U) 

The   local  truncation  error  in   eqs.    (l)    and   (2)    is   defined  as 

-y(tn+2}    +  yUn}    +   2h   y'(tn+l}    =   '    3  hV3)(tn)    +   0(hU) 
and 

-y(tn+2!  +  y(tn+1)  ♦  |(y<tn+1)  ♦  y'(tn+2»=  ^h3y<3'(tn)  ♦  o<h") 

respectively.      Hence,   if  y      and  y     _    are    correct,  we   expect  p    ,_  -  y 

n  n+1  n+2  n+2 

to  be    -  ^|-h3y^(t    )   +   0(h    )    as    in   eq.    (3).      Henrici    [3]   p.    257  also 

justifies  this   as   the  value   in  the   case   that   the   a.    and  a*  of  the 
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predictor   and  corrector  formulas   are   equal.      This  is   not  true   in  this 
example   and  we    find  that  eq.    (h)    does  not   agree  with  eq.    (3).      If  we 

subtract  eq.    (2)    from  eq.    (l),  we    find  that 

yp=y_y  +2hf(y  )_*Lf(y  )  (5) 

Jn+2        yn+2       'n       ^n+l        2        v,yn+ly        2      vjrn+2;  v?; 


In   this   paper  we    are   going  to  examine   the  nature  of   linear  combinations 
like  the   right-hand  side   of  eq.    (5)    under  the  assumption  that  y     is   the 
numerical  solution   derived  by   any  of  a  large   class   of  methods,    including 
multistep   and  one-step  methods. 


2.      HYPOTHESES 

We   assume  that   f(y)   has    continuous    (M+l)-st   derivatives   in   a 
suitable   region   so  that  the    initial  value  problem 

y'   =  f(y),  y(o)   given  (6) 

has    a  unique  solution  y(t)    for  t  e[0,T]  with  a   continuous  ■  (M+2)-nd 

derivative.      We   also   assume   that  we  have   a  technique   for   computing  a 

numerical  approximation  y     to  y(t    )   on   a  set   of  points 

n  n 

0=t0<   tx<   t2    ...    <   tH<T 

We  will  restrict    consideration  to  equally  spaced  points   so  that   t     =  nh. 

n 

However,    the   results    apply  to  the   unequal   spacing  given  by 

h     =  t   . _   -  t  (7) 

n  n+1  n 

with 

h  =     max     (h    )  (8) 

0<n<N       n 

when  the  hypotheses   of  this   section  hold. 

Many  numerical  techniques   are   such  that  the  numerical  approximation 

has   the   following  property    (see   Stetter    [5]). 

For  fixed  t 
n 

M 

y     =y(t    )   +     E      e.(t    )   hJ   +   0(hMfl)  (9) 

n  n  .  j      n 

J=P 

as  h  -*■  0,   n  ->  °° 
where  the  e.(t)   have   continuous    (M+l-j)-th  derivatives    and  satisfy 

J 

differential  equations   of  the  form 


el   =  f  e.   +   f  (10) 

j       y  j        j 

(in  eq.    (10),    as   in   all   future   equations,   explicit   dependence   of 

functions   on   t  is  not   shown   unless  necessary   for   clarity.      Unless 

otherwise   stated,    f  and  its    derivatives   are   evaluated  on  the 

solution  y(t).      Thus,    f     =    9f/9y|         ,,\.      The  subscript  n  will 

y  y=y ( t J 

refer  to  evaluation   at  the  point   (t    ,y    )  ,  while   the  superscript 

n     n 

n  will  refer  to  evaluation   at  the  point    (t    ,y(t    )).      Thus,   y     =  y(t    ).) 

n  n  n 

We  will  assume   that   eqs .    (9)    and  (10)   hold.      cj>      is    called  the 
principal  error  function   (see  Henri ci    [3],   particularly  Theorems    1.5> 

h.  3,    and  5. 12)    since   the   local  truncation  error  of  a  method  is 

p+1  /    P+2N  /    \ 

h  <{>     +   0(h        ).      In   general,    <t>.(t)    is   a  function   of  derivatives    of 

y   and  f  and  of  e.    for  i   <   j  ,   so  that   eq.    (10)    serves   to   define  the  e. 
i  J 

recursively.      In   a  multistep  method  in  which  sufficient    corrector 

iterations    are  used  or  the   order  of  the  predictor  is  high   enough, 

4>     =   C  y^P+      .      (See  Henrici    [3],   Section   5-3-7.) 

The  popular  Milne  method  for  estimating  local   errors    (Henrici    [3], 

Section   5-3-6)   has   been  seen  to  be   equivalent  to   forming  a  linear 

combination   of  values   of  y      .    and  computed  derivatives    f(y      .).      For 

n+i  n+i 

this  reason  we  wish  to  examine  the  linear  operator 


D,  [z]  =  E   [a.  z(t+ih)  +  b.hz' (t+ih)  ]  (ll) 

h      i=0  x 

(When  we  consider  unequal  steps,  we  must  consider 


D,  [z]   =  I   [a.  z(t  ..)  +  b.  hz'(t    )]  (lla) 

h   n    .  „   in   n+i     in     n+i 
i=0 

where   a.      and  b.      may  be   functions   of  h.) 
in  in 


For  sufficiently   differentiable   z,  we    can  write 

Djz]  =        E        d,z(j)hj   +  OCh^1)  (12) 

h  j=q+l      J 

(For  unequal  steps,  we    get 

M  /    » 

D.[z]     =        E        d.      z(JV   +  0(hMfl)  (12a) 

h        n        J-q+1      Jn 

where   d.      may   depend  on  h,   provided  that   a.      and  b.       are   sufficiently 
jn  in  in 

well  behaved  in   a  neighborhood  of  h   =  0.      This  we   must    assume.) 

In   fact,  when  we  estimate  the    error  of  a  numerical   solution,  we 

cannot   evaluate   eq.    (ll)    directly  as  we    cannot   assign   a  meaning  to 

y'.      Instead  we   must   evaluate   the   associated  nonlinear  operator 
n 

k 

D    [z]  =     Z      [a.    z(t+ih)   +  b.hf ( z+ih) ]  (13) 

i=0        X  X 

We   are   interested  in  the   asymptotic  expansion  of  D.  [z]  when    z  is 

identified  with  y      for   fixed  t    .      We  will   call  this   D,  [y    ].      Note  that 
n  n  h     n 

Dh[yn]    =  Dh[yn]  (1U) 

but   that  D,  [y    ]   is   not   defined, 
h     n 


3.      ASYMPTOTIC  ERROR  AND  DERIVATIVE  ESTIMATES 

In  this   section  we  will  expand  D    [y    ]   and  relate    it  to  D    [y    ]   and 

4>    .      This   result  will  show  how  y  or  4      can  be   estimated. 

P  P 

Substituting  eq.    (9)    into  eq.    (13)  we    get 

k  M  M 

Djy]=     E      [a.{yn+1  +      I      en+V }  +  b.hf  {yn+1  +     E      en+V}] 

h    n        i=o      1  j-p    J  1  3-p    J 

+  0(hMfl)  (15) 


We  write  the  Taylor's   expansion  of  f(y    )    about   f(y    )    as 

M 

f(yn)   =   f(yn+     E     eV)  +  0(hM+1) 

J=P 

=    fn  +      Z      fneV    +        E        rV   +   0(hMfl)  (l6) 

j=P     y   J  j=P+l     J 

where   r.    is   that   part  of  the   coefficient  of  h     which   involves   second 

and  higher   derivatives   of  f  with  respect  to  y.      Note  that   r.    depends 

also   on  e      for  k  <    j  .      Substituting  eq.    (l6)    into  eq.    (15)    and  defining 

K. 

r     =  0  we   get 
P 

W  -  \  [(vn+i  +  v^1' 

1=0 

M 

+      E      (a.en+i  +  b.h{fn+i   en+i  +   rn+i})hJ]  +   0(hMfl)  (lT) 

J=p        1   J  x        y  j  j 

We   substitute   f  e.    =  e!    -   d> .    from  eq.    (10)    into  this   to   get 
y   J  J  J 

M 

L\[y    ]   =   D    ryn]   +      Z      D Jen]hJ 
h     n  h  .  h      j 

J=P 

M-i    k  . 

-     E        E      b.?fV+1  +   0(hMfl)  (18) 

j=P  i=0      X    J 


where 


> .    =    <p .    -   r . 
J  J  J 


(19) 


The    continuity  hypotheses    allow  h  <f> .        to  be  expanded  by  Taylor's 

J 

series   about   t      as 
n 

M-i-1  (t    ) 

hJ+1  ?+i  =    I    hJ+l+1  J*  JLsl  ♦  oth*1)                     (20- 

0               t=0                        J  I. 


where 

t   .  -  t 

T  =  -=-= —  (=  i  for  equal  intervals) 

i       h 

When  eq.  (20)  is  substituted  into  eq.  (l8)  we  get  the  desired  result: 

M 

W  -  V^>  +  ,E  Dh[ej!  h° 

-  V  hJ+1  ¥  B,    ?(iJn  ♦  0(hW1)  (21) 

where 

k     b.(T.)£ 

\  ■    =    -TT-  (22) 

i=0 
From  eq.    (12)    and  the   continuity   assumptions, 

D.[en]hJ   =   Odi1*11^1'    J+q+l))  (23) 

Eqs.  (21)  and  (22)  allow  us  to  make  the  following  observations: 
Observation  1 

If 

(i)    M  >  q  <  p 

or        (ii)   M  >  q  >  p  and  B_  =  B   =...=  B     =0 

—        0     1        q-p 

then  Dh[yn]   =   dq+1  h^1  y(q+l)n  +  0(h<1+2)  (2k) 

Observation   2 


If  q  =  p   <   M,    then 

f>h[yn]  =  (d        y(<1+1)n  -  BQ  fy**-*1  +  o(h*+2)  (25) 


Observation  3 


If  q  >   p   <   M  and  B     4  0 ,   then 

W  -  -\  *p  *p+1  +  °<»p+2> 


k.      APPLICATIONS 

The   result   in  the  previous    section  will  be  applied  to  justify 
two   common   techniques    and  one  new  one  for   error  estimation   in  this 
section.      Note  that   the  numerical  solution   is   still  assumed  to   satisfy  eq.    (9) 
Lemma  1 

Differences    can  be   applied  either  to  the  solution  y     or  to 
the    computed  derivatives    f(y    )    to   estimate  y  for   q   <   N. 

Proof 

This    follows    directly   from  observation   1.      If  the    differences 

of  y      are   formed,   then   all  b.    and  hence  Bn     are   zero.      If  differences 
n  1  £ 

of   f(y    )    are   formed,   it   is   trivial  to  verify  that  B     =  0   for 
n  >v 

0    <_  l     '-_  q-l.      Since  p   >_  1   for  eq.    (9)    to  be   meaningful,    one   of  the 

alternatives   of  observation   1   applies. 

Q.E.D. 

This   means,    for  example,   that  the   common   technique   of 

estimating  the   error  in  Adam's    method  is      valid  since   backward 

differences   of  f(y    )    are   computed.      These   can  be   used  to   estimate 

the  error  in   any  order  method,   not    just   in   the   order  used  to   compute 

the  y    .      A  simple  extension   can   also  be   used  to   justify  the   following 

technique   for   fixed  steps: 

Suppose  we  have   an  estimator  D,  [y    ]    for  h  y  which 

h     n 

satisfies   the   conditions   of  observation   1.      Then  Dn  [y    ]   -   D.  [y     .  ] 

h     n  h     n-1 

q+2       (q+2) 
is    an   estimator   for  h  y  ifq<    M-l.      This    follows   by 

considering  the   operator  D    [•     ]   -  Dn  [•       ,  ]  =   fi,  [y    ].      If  the    coefficients 

hn  hn-1  hn 

B„     correspond  to   f)     and  B„     to  Dn  ,    then  Bn  =   0.      If  B .    =   0   for 
I  h  I  h  0  j 

0   <_  j  £  q-p,   then  B.    =  0   for  1  <_  j   <_  q-p+1 .      This  technique   is 

J 

frequently  used  to   estimate   the   error  of  a  multistep  method  of  order 
one   larger  than  that   currently   in  use.      D    [y    ]   is   the  predictor- 
corrector  difference    (see   Lemma  2  below)    and  it    is  differenced  again. 
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Lemma  2 

If  a  P(EC)   v  predictor-corrector  method  is  used  and  the  predictor 
and  corrector  have   the   same   order  p,   then   the   difference  between   the 
predicted  and  corrected  value   is   a  multiple   of  h       '  y  P  +  0(h        ) . 

Proof 

If  the  predictor  is 

pn+k=   .E      <«i  Vi  *  h»i   f(yn+i"  <27) 

1=0 

and  the   corrector  is 

U  =    .\    (ai   Vi   +  hBi    f{W   +  hBk   f(W  (28» 

1=0 

then  D,  [y   ]  =  p         -  y 

h     n  n+k  n+k 


=      E      a.    y  +  hb.    f(y        ) 

.    _      i     n+i  l  n+i 

i=0 

where   a.    =  a.    -  a*  and  b.    =   3.    -3*   (assuming  3n     =  0   and  a,     =  a  *  =   -l) . 
111111  k  kk 

If  the   local   errors    in  the  predictor  and  corrector  are   defined  as 

E      [a.    y(t    ..)    +  h3.    y'(t    ..)]  =  c    ._    hP+1  y(P+l)n   +  0(hP+2)        (29) 
.    ~        i  n+i  l  n+i  p+1 

i=0 


and 

k 

E 

i=0 

then 


E      [a*  y(tn+.)   +  h3*  7f(tn+1)]   =  c*+±  hP+1  y(p+l)n  +  0(hP+2; 


(30) 


Dhty]   =    (cp+1  -   c*+1)   hP+1  y(p+l)    +  0(hP+2)  (3D 

In  this    case    (see  Henrici    [3],   eq.    (5-185)) 
JB*L 


*      =     k  (32) 

p         E      3* 
i=0     X 

so   from  observation  2, 
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k-1 

£      6. 
i=0      X 


p  y  =     fi      [y      ]     =      [C  ~ C*  ]     hP+1     y(^)»     +     0  (h^  )       (33) 

^n+k  n+k  h      n  p+1  ^     p+1 

i=0     i 

Q.E.D. 

Comment 

This   result   is   an  extension  of  the  Henrici    result   mentioned  in 

Section   1,    and  also   indicates    the   correction   factor  to  be   used  when 

k-1  k 


Z     B.    ¥     £     3*. 
i=0     x       i=0     x 

Example 

For  the  example 

in 

c3  -  -i,         Mi 

=  2 

c*  =  i.         ZBf 

=  1 

The   Milne  estimate   suggests   that  the  predictor-corrector 
difference  should  be 

(c3  -  c*)   h3  y(3)   +  0(hU)    =  ^h3  et  +   0(hU) 

which   is   true    for  p     -  y     only.      But  this   is    for  n   fixed,  not    for  t 
fixed.      In  the  latter   case,   eq.    (33)    tells   us   that 

*n+2->V2   =  -Jh3^(3)   +   0(hH)  <*> 

This   example  points   out   another  problem.      Why,    one   should  ask, 
should  we  not  set  n  =  0   in  eq.    (3*0    to   get 

P2   -  y2  =   -  |h3  +  0(hS 

which  has  been  seen   to  be   incorrect,      t      is,    after   all,    a  fixed  value. 

The   answer  is  that  the   initial  hypothesis   implied  by  eq.    (9)    is  not 

valid  at   t     =0    for  any   M  >   1,  but   it   is   valid  at   t      >   0   for    aiy  M. 
n  n 
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This   is   due   to  the  starting  "errors"  which    arise  when   additional  values 
are  used  to  start  multistep  methods.      The  corrector  eq.    (2)    is   a 
special   case   of  a  multistep  method — namely   a  one-step  method.      In 
fact,   the   corrector  is   started  at   t     =  h.      It   is   simple   to  verify  that 
for  t   >   0,   e     =  —  te      and  e_  =  -  rr  e    (2-t)    if  this    corrector  is 
solved  exactly.    However,   eq.    (9)    does  not  hold  at  t  =  0  with   these 
values   since  e    (0)    f  0.      Viewed  as   a  two-step  method,   eq.    (U)   has 
the   associated  stability  polynomial  p(0    =   -£     +    £.      Its   nonprincipal 
root. is   zero.      In   the   case   of  a  general  multistep  method  with  nonzero 
nonprincipal  roots,   starting  errors   may   continue   to  have   an   effect. 
Henrici    [3],   Section   5.3-5   discusses   this    and  shows  that  the  starting 
errors   do  not   decay  for  weakly   stable  methods.       (See   eq.    (5-205), 
ibid.)      Nonprincipal  roots    £.    introduce   terms    of  magnitude    (£.)    e 

into  y     -  y     where   e   is   the   magnitude   of  the  starting  error.      For 

i       i  /    M+-l\ 

I  £.  |    <   1,   these   are   smaller  than  0(h        J    for  any  M,  but   remain  of 

0(e)    for    |  E,.  |    =1.      Consequently,  we   cannot  expect  weakly   stable 

methods   to   satisfy  eq.    (9)    unless  we  have    a  fortunate    (and  unlikely) 

choice  of  starting  values.      Strongly   stable  methods  will  satisfy 

eq.    (9)    for  any  fixed  t      away      from  the  starting  region.      Thus,   in 

the   example   above,    "e    (t    )"   is   actually   given  by 

*3(V  " "  &  •*"  <2-V  +  k  (52>n  .     (35) 

for   £     =  0.      This   fails  to   satisfy  the   differentiability  properties 
at  t  =  n   =  0 . 

This   problem  has   serious   implications   since   it   means   that  the 
estimate   of  the  error  immediately  following  a  change   in   step   size   or 
order  may  not  be  valid. 
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Our  final  application   is   to  methods   in  which  the   local  error 

d>     h  is    a  complicated,  combination   of  derivatives   which   cannot 

P 

reasonably  be    computed.      In   this    case,  we   use   observation    3.      If  a 

differentiation   formula  EL  [z]   is    chosen   for  which  B     ^   0,    and  for 

which   q  >   p,  we   get   a   direct   estimate    of  B„   4    . 
^        r  Op 

Example 

Consider  the  Heun  method 

k0  "   hf(yn» 

kl  "  hf(^„  +  V 

This   is   a  two-stage  Runge-Kutta  method.      The   asymptotic  local 
truncation   error  is 

h3  *2  =  h3(y(3)/12  -   fy  y(2)A)  (37) 

The  numerical  solution   satisfies   eq.    (9).      If  we    consider 

D    [z]   =  -z(t+2h)    -   Mt+h)   +   5z(t)   +   Uhz'(t+h)   +  2hz'(t) 

--$hVfc)    +0(h5) 

we  expect   to    find  that 

Vyn]    =  "B0    *2   h3  =    "6    ^2   h3'  (38) 

Consider  y'    =  y,   y(0)    =   1.      It   follows   immediately  that 

,2 

y;  =  yn  =  (i  +  h  +  |r)n  (39) 

and  that 

*2  =  -y/6  Uo) 

If  we   evaluate   D.  [y    ]    for  eq.    (39)  ,  we    get 
h     n 


Ik 


Dh[yn]  =  h3  yn  +   0(hU) 

agreeing  with  eqs.    (38)    and   (Uo). 

Example 

The  error  in   any   fourth  order  Runge-Kutta  method  could  be 
estimated  using 

D    [z]  =   -z(t+3h)    -  l8z(t+2h)   +  9z(t+h)   +   10z(t)   +   9hz(t+2h) 

+   l8hz(t+h)   +   3hz(t)   =   0(h)  (Ul) 

It    gives 

*U  =   "Dh[yn]/3°  +  °(h6) 

In  practice,   this   requires   values    at   two  preceding  points   to 
estimate  the  error  in  the  step  just  performed.      This   probably   limits 
the   flexibility  of  eq.    (Ul)    as   an  estimator  since   those    steps   have 
to  be  of  constant   size.      The  equivalent  of  eq.    ( Ul )    could  be   obtained 
for  variable   steps,   but  eq.    (9)    does   not  hold  under  most    circumstances, 
so  the  estimate  may  not   be  valid. 

5.      CONCLUDING  REMARKS 

Little  has   been   said  about  the    conditions  under  which  eq.    (9)    holds, 
as    a  survey  of  all  types   of  methods  would  occupy  a   considerable    amount 
of  space.      We  will  summarize   the  situation  briefly   and  indicate   some 
open  areas. 

In  the   case   of  constant  steps,   eq.    (9)    holds    for  any  one   step 
method  whose   local  truncation   error  has    an   asymptotic  expansion.      (This 
includes   all   one  step  methods  known   to  the  writer.)      Eq.    (9)    is    also 
valid  for  strongly  stable  P(EC)n±l   and  corrector  only  multistep  methods 
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(but  note  the   starting  problems    indicated,  by  example   in   the  last 

section).      P(EC)      methods   present  problems  which  need  to  be   examined 

because   there   are    two    different  values    of  y      saved  from  step  to   step — 

the   final  corrected  value  y     and  the  penultimate  value  y     used  in  the 

n  n 

last   evaluation   of   f(y    ). 

n 

Variable  steps  present  serious  problems.   If  there  exists 

a  sufficiently  smooth  function  9  (t)  such  that 

0  <  A  <_  e(t)  <  1  (1+2) 

and  h  =  hO(t  ) 
n       n 

it   is  possible  to  show  that   the  above  methods   satisfy  eq.    (9)-      For 

one  step  methods  this   is   an  extension  of  the  approach  used  in 

Henrici    [3]    and  Gear    [l],   Section    U.7.      For  multistep  methods   it   is 

necessary  to  verify  that  Stetter's    [5]   results   can  be   applied.      Gear 

and  Tu   [2]   show  that  eq.    (U2)    is   sufficient   for  stability   for 

"reasonable"   variable   step  multistep  methods.      These    are   multistep 

methods  whose    coefficients    can  be   expressed  via  eq.    (U2)    as   functions 

of  h  which  are   sufficiently   smooth  in   a  neighborhood  of  the  origin. 

The   analysis    seems  to  have   little  meaning  unless   it   is   possible 

to  express   all   step   sizes  h      as   smooth  functions   of  h,    and  it   is 

questionable  what  this   means   in   the   context   of  a  computer  program! 
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